/*		
	This header-file is filled with functions from solutions to previous assignments:
	               (behind '~' the motivation for including these functions is stated)
	-graymorphology.h		~ for erision and dilation
	-MoreFiltering.h		~ for canny-edge / gaussian blurring
	-filter.h & filter.cpp	~ for applying canny-edge / gaussian blurring
	-rgh.h                  ~ for pixel-class
	The functions were created by supervisors of the image processing courses.
	As they are exemplarely solutions, only changed lines of code are commented.
	Some unneeded functionalities were removed without comment.
	The un-templated function bodies are to be found in 
	-ip.cpp.
*/

#ifndef __IP_H__
#define __IP_H__

#include <algorithm>
#include <iostream>
#include <string>
#include <vector>
#include <queue>
#include <list>
#include <cmath>

#include "Image.H"
#include "IO.H"
#include "Util.H"

/**
 * Returns the value at (x, y) using linear interpolation. See also
 * Sonke, page 122 and 123, equation 5.22
 *
 * @param image Input image
 * @param x
 * @param y
 * @return Interpolated value at (x, y)
 */
float linearInterpolation(const pip::Image<float> &image, const float &x, const float &y);
/**
 * Convolves the input image with a kernel and stores the result in the output.
 *
 * @param input Constant reference to the input image
 * @param output Reference to the output image
 * @param mask Constant reference to the kernel
 */
void filter(const pip::Image<float> &input, pip::Image<float> &output, const pip::Image<float> &mask, const pip::BorderInfo &bi);
/**
 * Performs a separable convolution. The first dimension is convolved with the first
 * kernel in kernels, and so on.
 *
 * @param in Input image
 * @param out Output image
 * @param kernels List of kernels
 * @param bi Border information
 */
void filterSeparable(const pip::Image<float> &in, pip::Image<float> &out, const std::vector<std::vector<float> > &kernels, const pip::BorderInfo &bi  = pip::border_mirror);
/**
 * Samples the unnormalized Gaussian for given sigma at integer values.
 *
 * @param i
 * @param sigma
 * @return Unnormalized Gaussian
 */
float g0(const int &i, const float &sigma);
/**
 * Samples the unnormalized first derivative of the Gaussian for given sigma at integer values.
 *
 * @param i
 * @param sigma
 * @return Unnormalized first derivative of the Gaussian
 */
float g1(const int &i, const float &sigma);
/**
 * Samples the unnormalized second derivative of the Gaussian for given sigma at integer values.
 *
 * @param i
 * @param sigma
 * @return Unnormalized second derivative of the Gaussian
 */
float g2(const int &i, const float &sigma);
/**
 * Samples the unnormalized third derivative of the Gaussian for given sigma at integer values.
 *
 * @param i
 * @param sigma
 * @return Unnormalized third derivative of the Gaussian
 */
float g3(const int &i, const float &sigma);
/**
 * Normalizes the vector to unit absolute sum.
 *
 * @param kernel Kernel to normalize
 */
void normalize(std::vector<float> &kernel);
/**
 * Computes a Gaussian derivative of scale sigma and order dx, dy of in, result stored in out.
 *
 * @param in Input image
 * @param out Output image
 * @param sigma
 * @param dx
 * @param dy
 * @param cutoff
 * @param bi Border information
 */
void gauss(const pip::Image<float> &in, pip::Image<float> &out, const float &sigma, const int &dx = 0, const int &dy = 0, const float &cutoff = 5.0 , const pip::BorderInfo &bi = pip::border_mirror);



/// CANNY EDGE DETECTION \\\

/**
 * Gradient magnitude ...
 *
 * @param in Input image
 * @param dx Output of X-derivative
 * @param dy Output of Y-derivative
 * @param gm Output of Gradient Magnitude
 * @param sigma
 */
template <class T> // also this function is now templated...
void gradientMagnitude(const pip::Image<T> &in, pip::Image<float> &dx, pip::Image<float> &dy, pip::Image<float> &gm, const float &sigma)
{
  Image<float> image;
  gm.resize(in.dimension());

  convertAllToFloat(in, image);

  gauss(image, dx, sigma, 1);
  gauss(image, dy, sigma, 0, 1);

  for(int y = 0; y < gm.dim(1); ++y) {
    for(int x = 0; x < gm.dim(0); ++x) {
      gm(x, y) = sqrt(dx(x, y) * dx(x, y) + dy(x, y) * dy(x, y));
    }
  }
}

/**
 * Canny edge detector ...
 *
 * @param in Input image
 * @param out Output of gradient magnitude for pixels that survive non-maximal suppression
 * @param sigma Scale for the gradient magnitude
 */
template <class T> // adding template to operate on shot images instead of unsigned char only
void cannyEdgeDetector(const pip::Image<T> &in, pip::Image<float> &out, const float &sigma)
{
  out.resize(in.dimension());
  out.fill(0);

  pip::Image<float> dxImage, dyImage;
  pip::Image<float> gmImage(in.dimension());

  // Compute gradient magnitude and edge direction
  gradientMagnitude(in, dxImage, dyImage, gmImage, sigma);

  // Apply non-maximum suppression in direction of the edge
  for(int y = 0; y < in.dim(1); ++y) {
    for(int x = 0; x < in.dim(0); ++x) {
      float gMag(gmImage(x, y));

      //if(gMag > 1.0E-10) { //if smaller, set to 0
        // calculate normal to edge according to formula 4.63 in Sonka
        // normal = derivative/gradientmagnitude
        float dx(dxImage(x, y) / gMag),
          dy(dyImage(x, y) / gMag);

        // do non-maximal suppression --> get only the strongest edge
        // Get gradientmagnitude interpolation value in positive and negative direction normal to edge
        // to check if these edges are stronger!
        float gMag1(linearInterpolation(gmImage, x + dx, y + dy)),
          gMag2(linearInterpolation(gmImage, x - dx, y - dy));

        if(gMag1 < gMag && gMag2 < gMag) { // x, y is maximum
          out(x, y) = gMag;
       // }
      } else {
        out(x, y) = 0;
      }
    }
  }
}

//
//   code of graymorphology.h:
//
template <class T> // templating this function to operate not only on char images
void discStructureElement(pip::Image<T> &imOut, int rad){
     pip::Dimension dim(2*rad+1,2*rad+1); // dimensions of SE --> +1 for center
     imOut.resize(dim);
     int rad2 = rad*rad;
     int val;
     for (int y=0; y<=rad; ++y){
        for (int x=0; x<=rad; ++x){
          val = 0;

          if( (x*x+y*y) < rad2+0.5)
              val = 1;

          imOut(rad+x,rad+y) = T(val);
          imOut(rad-x,rad+y) = T(val);
          imOut(rad-x,rad-y) = T(val);
          imOut(rad+x,rad-y) = T(val);
      }
	}
}

template <class Tin>
void grayDilate(const pip::Image<Tin> &imIn, pip::Image<Tin> &imOut, const pip::Image<Tin> &se){
        // make sure we have the correct output size
        if (imOut.dimension()!= imIn.dimension()) {
                imOut.resize(imIn.dimension());
        }

        int iSEWidth  = se.dim(0);
        int iSECenter = iSEWidth/2;
        int iWidth    = imIn.dim(0);
        int iHeight   = imIn.dim(1);

        int iX, iY;
        // over all image pixels
        for (int y=0; y<iHeight; ++y)
        for (int x=0; x<iWidth; ++x) {
                imOut(x,y) = imIn(x,y);
                // over all se pixels
                for (int m=0; m<iSEWidth; ++m)
                for (int n=0; n<iSEWidth; ++n) {
                        if (se(n,m) > 0) {
                                iX = x+(m-iSECenter);
                                iY = y+(n-iSECenter);
                                if ( iX >= 0 && iX < iWidth && iY >= 0 && iY < iHeight)
                                imOut(x, y) = max(imOut(x,y),imIn(iX,iY));
                        }
                }
        }
}

template <class Tin>
void grayErode(const pip::Image<Tin> &imIn, pip::Image<Tin> &imOut, const pip::Image<Tin> &se){
        // make sure we have the correct output size
        if (imOut.dimension()!= imIn.dimension()) {
                imOut.resize(imIn.dimension());
        }

        int iSEWidth  = se.dim(0);
        int iSECenter = iSEWidth/2;
        int iWidth    = imIn.dim(0);
        int iHeight   = imIn.dim(1);

        int iX, iY;
        // over all image pixels
        for (int y=0; y<iHeight; ++y)
        for (int x=0; x<iWidth; ++x) {
                imOut(x,y) = imIn(x,y);
                // over all the pixels
                for (int m=0; m<iSEWidth; ++m)
                for (int n=0; n<iSEWidth; ++n) {
                        if (se(n,m) > 0) {
                                iX = x+(m-iSECenter);
                                iY = y+(n-iSECenter);
                                if ( iX >= 0 && iX < iWidth && iY >= 0 && iY < iHeight)
                                imOut(x, y) = min(imOut(x,y),imIn(iX,iY));
                        }
                }
        }
}


// Return minimum and maximum of the input image
template <class T>
void getMinMaxImage(const pip::Image<T> &in, T &min, T &max)
{
  min = numeric_limits<T>::max();
  max = numeric_limits<T>::min();

  for (int z = 0; z < in.dim(2); ++z) {
    for (int y = 0; y < in.dim(1); ++y) {
      for (int x = 0; x < in.dim(0); ++x) {
        T value(in(x, y, z));

        if (value > max) {
          max = value;
        }

        if (value < min) {
          min = value;
        }
      }
    }
  }
}

// Arrays with offsets to 2D neighborhood.
// Taking the first 4 gives 4-neighborhood, taking all gives 8-neighborhood.
static const int dnbx2D[8] = {1, 0,-1, 0, 1, 1,-1,-1};
static const int dnby2D[8] = {0, 1, 0,-1, 1,-1, 1,-1};
// Arrays with offsets to 3D neighborhood.
// Taking the first 6 gives 6-neighborhood, taking first 18 gives 18-neighborhood,
// taking all gives 26 neighborhood.
static const int dnbx3D[26] = {1, 0, 0,-1, 0, 0,   1, 1, 1, 1, 0, 0,-1,-1,-1,-1, 0, 0,   1, 1, 1,-1,-1,-1, 1,-1};
static const int dnby3D[26] = {0, 1, 0, 0,-1, 0,   1,-1, 0, 0, 1, 1, 1,-1, 0, 0,-1,-1,   1,-1, 1, 1, 1,-1,-1,-1};
static const int dnbz3D[26] = {0, 0, 1, 0, 0,-1,   0, 0, 1,-1, 1,-1, 0, 0, 1,-1, 1,-1,   1, 1,-1, 1,-1, 1,-1,-1};



class Pixel
{
        private:
            int x, y, z;

        public:
        Pixel(int x = 0, int y = 0, int z = 0) {set(x,y,z);};

        template <class T>
        T get(const pip::Image<T> &im) const {return im(x,y,z);}

        template <class T>
        void set(pip::Image<T> &im,T val) const {im(x,y,z) = val;}

        void get(int &ix, int &iy, int &iz) const {ix=x;iy=y;iz=z;}

		bool negative()  {return (x < 0);}                            // added function to identify marker-pixel (out of image-range) 

        void set(int ix, int iy, int iz = 0) {x=ix;y=iy;z=iz;}

		Pixel& operator=(const Pixel &rhs) {if (this != &rhs) { x=rhs.x; y=rhs.y; z=rhs.z; } return *this; }

        template <class T>
		void getNeighbors(const pip::Image<T> &im,std::vector<Pixel> &v, int n=6) const
        {
                v.clear();
                if(n==4||n==8){//2D
                        for(int i=0; i<n; ++i){
                                Pixel p = *this;
                                p.x+=dnbx2D[i];
                                if(p.x<0 || p.x>=im.dim(0)) continue;
                                p.y+=dnby2D[i];
                                if(p.y<0 || p.y>=im.dim(1)) continue;
                                //if we get here, the position is in the image
                                v.push_back(p);
                        }
                }
                else if(n==6||n==18||n==26){//3D
                        for(int i=0; i<n;++i){
                                Pixel p =*this;
                                p.x+=dnbx3D[i];
                                if(p.x<0 || p.x>=im.dim(0)) continue;
                                p.y+=dnby3D[i];
                                if(p.y<0 || p.y>=im.dim(1)) continue;
                                p.z+=dnbz3D[i];
                                if(p.z<0 || p.z>=im.dim(2)) continue;
                                v.push_back(p);
                        }
                }
                else
                assert(false);//illegal value of n
        }

};


#endif